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We apply adjoint-based sensitivity analysis to a time-delayed thermo-acoustic system: 
a Rijke tube containing a hot wire. We calculate how the growth rate and frequency 
of small oscillations about a base state are affected either by a generic passive control 
element in the system (the structural sensitivity analysis) or by a generic change to its 
base state (the base-state sensitivity analysis) . We illustrate the structural sensitivity by 
calculating the effect of a second hot wire with a small heat release parameter. In a single 
calculation, this shows how the second hot wire changes the growth rate and frequency 
of the small oscillations, as a function of its position in the tube. We then examine the 
components of the structural sensitivity in order to determine the passive control mecha- 
nism that has the strongest influence on the growth rate. We find that a force applied to 
the acoustic momentum equation in the opposite direction to the instantaneous velocity 
is the most stabilizing feedback mechanism. We also find that its effect is maximized 
when it is placed at the downstream end of the tube. This feedback mechanism could be 
supplied, for example, by an adiabatic mesh. We illustrate the base-state sensitivity by 
calculating the effects of small variations in the damping factor, the heat-release time- 
delay coefficient, the heat-release parameter, and the hot wire location. The successful 
application of sensitivity analysis to thermo-acoustics opens up new possibilities for the 
passive control of thermo-acoustic oscillations by providing gradient information that can 
be combined with constrained optimization algorithms in order to reduce linear growth 
rates. 



1. Introduction 

In a thermo-acoustic system, heat release oscillations couple with acoustic pressure 
oscillations. If the heat release is sufficiently in phase with the pressure, these oscillations 
grow, sometimes with catastrophic consequences. Using adjoint sensitivity analysis, we 
identify the most influential components of a thermo-acoustic system and quantify their 
influence on the frequency and growth rate of oscillations. This technique shows how a 
thermo-acoustic system should be changed in order to extend its linear ly stable reg ion. 



Adjoint sensitivit y analysis of incompre s sible flows was proposed by iHill (1992) and 



developed further bv lGiannetti fc Luchini I ( 2007 ) in order to reveal the region of the flow 



that causes a von Karman vortex street behind a cylinder. They used adjoint methods 
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to calculate the effect that a small control cylinder has on the growth rate of oscilla- 
tions, as a function of the control cylinder's position downstream of the main cylin- 
der. This control cyli n der in duce s a force in the opposite direction to t he velocity field. 
Giannetti fc Luchini (j2007t) and iGiannetti. Camarri fc Luchini J 2010l) c onsidered this 



feedback only on the perturbed fields but lMarquet. Sipp &: Jacauin J 20081) extended this 



analys is to consider the cylinder's effect on the base flow as well. lSipp. Marquet. Meliga fc Barbagallo 



(2010) p rovide a comprehensive review of sens itivity analysis for incompressible fluids and 
IChandler. Juniper. Nichols fc Schmid I (|2012t ) extend this analysis to low Mach number 
flows in order to model variable density fluids and flames. 

The aim of this paper is to extend adjoint sensitivity analysis to a thermo-acoustic 
system, whic h has not been attempted before. We inve stigate the therm o-acoustic system 
described by iBalasubramanian fc Suiith I ( 2008al) and Juniper ( 2011 ). This is an open- 
ended tube, through which air passes, and which contains a hot wire at a given axial 
location. One-dimensional acoustic standing waves in the tube modulate the air velocity 
at the wire, which in turn modulates the heat transfer from the wire to the air, which 
is modelled with a modified form of King's law (|Heckl 19901 Matveev 2003). This heat 
transfer occurs at the wire's location but is not instantaneous. The time taken for the heat 
to diffuse to the bulk fluid is modelled as a time delay between the velocity fluctuations 
and the heat release fluctuations. 

The analysis consists of three main steps. Firstly, we study the system as an eigenvalue 
problem in the complex frequency domain. Secondly, we derive two sets of adjoint equa- 
tions from the linearized governing equations. Thirdly, we use the adjoint equations to 
perform both a structural sensitivity analysis and a base-state sensitivity analysis. The 
structural sensitivity analysis quantifies the effect that feedback mechanisms have on the 
frequency and growth rate of oscillations. This analysis relies on studying the effect of 
a perturbation to the governing equations, which is known as a structural perturbation. 
There are several components of the feedback and, in this paper, we calculate all of them. 
We then illustrate the structural sensitivity by considering the effect of feedback from a 
second hot wire. The base state sensitivity analysis quantifies the effect of a change in the 
constant coefficients of the governing equations. It does not involve a feedback mecha- 
nism. The base state in this thermo-acoustic model is represented by four parameters: the 
damping factor, the heat-release time-delay coefficient, t; the heat-release parameter, 
f3, and the hot wire location, xh- This shows us how to change these parameters in order 
to most stabilize the system. In addition, we can also calculate the location of the first 
hot wire that makes the system most sensitive to base-state modifications. In the final 
section we apply this analysis to the passive control of an unstable nonlinear system. 



2. Thermo-acoustic model 



The thermo-acoustic system examined in this paper is a horizontal Rijke tube contain- 
ing a hot wire. It is governed by the following nonlinear time-delayed equations: 



= 0, 




+ u(t - t) 



S(x - x h ). 



(2.1) 
(2.2) 

(2.3) 
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where u, p and q are the non-dimensional velocity, pressure, and heat-release rate, re- 
spectively. The hot wire is placed at i = xh, which is modelled by the Dirac delta 
(generalized) function S(x — Xh). The system has four control parameters: £, which is 
the damping; /3, which encapsulates all relevant information about the hot wire, base 
velocity, and ambient conditions; r, which is the time delay, and Xh, which is the position 
of the hot wire. The values of (3, t, and Xh are given in the figure captions along with 
the damping constants c j and eg . In we will explain how C, is related to c\ and c 2 . 
Eqs. (|2.1[) - (l2.2j) are derived from th e Navier-Stoke s and energy equations by assuming 
first-order acoustics, as explained in ICulick I (119711). The heat-release rate in eq. (I2.3P is 
modelled with a modified form of King's law (|Heckl1ll990t iMatveev I [20031) . Note that 
throughout this paper we defin e the hea t -relea se parameter j3 to be v3/2 times the heat- 
release parameter (3 defined in lJuniper I (|201ll) . The heat- release term (|2.3|) is linearized 
around a fixed point of the system, where \uh\ <C 1. In addition, eq. (|2.3I) is linearized 
also in time assuming that the time-delay coefficient is sufficiently small compared with 
the period of the highest Galerkin mode (Sj3|: 



q = I ( U - ) (V(.r .17, i. 



(2.4) 



By substituting eq. (|2.ip into eq. (|2.4I) . we obtain an equivalent expression for the lin- 
earized heat-release law: 



P 



Xh)- 



(2.5) 



It is important to anticipate that, although eq. (12.41) is physically equivalent to (I2.5[) . the 
systems of the linearized governing equations ([2TT|) - ([2~2l) - (l2T4]) and ([2"TT ]I - ([%% ]) - ([2"3 |) will 
produce two different sets of adjoint equations (@. 



3. Numerical discretization 

The partial differential equations (|2.ip - (|2.2[) - (l2.4[) . which govern the thermo-acoustic 
system, are discretized into a set of ordinary differential equations by choosing an or- 
thogonal basis that matches the boundary conditions. This procedure is also known as 
the Galerkin method. The variables are expressed as: 



N N 



u(x,t) = 2_^r) j (t)cos(jirx), p(x, t) = - ^ I — — 1 sin(jTnc). (3.1) 
i=i j=i ^ J7T ' 

The state of the system is given by the amplitudes of the Galerkin modes that rep- 
resent velocity, T)j, and those that represent pressure, rjj/jn. The state vector of the 
discretized system is the column vector x = ( u iP) 7 \ where u = (rji , . . . ,?yjv) T and 
p = (jii /n, . . . ,rjff I Nn) T . The discretized problem can be represented in matrix nota- 
tion: 

f - rx- 0.2) 

where f is the 2Nx2N direct matrix and x is the 2N xl state vector. The basis functions, 
cos(j7ra;) and sin(j7rx), are the eigenfunctions of the undamped acoustic system without 
the heater. The direct matrix F is shown in appendix \K\ in eq. (|A 1[) . Note that, when 
the system has N Galerkin modes, it has 2N degrees of freedom. 

The linearized equations in [J2] are valid for small \v,h\ and r <C T^, wh ere Tj = 2/j 
is the period of the j th Galerkin mode, as explained in lJuniper" ( 2011 ). T he results 



are presented here for a system with 10 Galerkin modes (as for system C in Juniper I 



4 L. Magri, M.P. Juniper 

(2011)). We checked modal convergence considering more Galcrkin modes and found 
that 10 modes provide an accurate representation of the system, as discussed in ^ 3T.31 

At the ends of the tube, p and du/dx are both set to zero, which means that the system 
cannot dissipate acoustic energy by doing work on the surroundings. Dissipation and end 
losses are modelled by the damping parameter for each mode Q = c\j 2 + C2\/~j, where 
ci and C2 are constants. Oscillations of higher Galerkin mo des decay very rapidly if no 



mechan ism drives them. This damping model was us ed in Balasu bramanian fc Suiith 
(2008b) and was b a sed o n correlations developed by iMatveev I (|2003h from models in 



Landau fc Lifshitz I (Il959h . 



4. Adjoint operator 

In this section the adjoint operator i s defined. This definition is an extension over the 
time domain of the definition given by Dennerv fc Krzvwicki ( 199(3). Let L be a partial 



differential operator of order M acting on the function q(xi , X2, ■ ■ ■ , xk, t), where K is 
the space dimension, such that Lq(xi ,x%,...,XK,t) — 0. We refer to the operator L as 
the direct operator and the function q as the direct variable. The adjoint operator L + 
and adjoint variable q + are defined via the generalized Green's identity: 



J J <tLq- q (L+q+)dVdt = J J J2 nidSdt+J Qi(q,t)\odV. 

(4.1) 

where i = 1,2, . . . , K and Qi(q,q + ) are functions which depend bilinearly on q, q + 
and their first M — 1 derivatives. The complex-conjugate operation is labelled by an 
overline. The domain V , is enclosed by the surface S, for which rij are the projections on 
the coordinate axis of the unit vector in the direction of the outward normal to the surface 
dS. The time interval is T. The adjoint boundary conditions and initial conditions on 
the function q + are defined as those that make the RHS in eq. (|4.1|) vanish identically on 
S, t = and t = T. 

The adjoint equations can either be derived from the continuous direct equations and 
then discretized (CA, discretization of the Continuous Adjoint) or be derived directly 
from the discretized direct equations (DA, Discrete Adjoint). For the CA method ( ^6.11 
and fc|6.2p . the adjoint equations are derived by integrating the continuous direct equations 
by parts and then applying Green's ident ity (14.11). They are then discretized with the 



Galerkin method (|3.1[) . The appendices of lJuniper I ([20111 ) show the intermediate steps 



Two different sets of adjoint equations are derived here, shown in table [TJ The first 
set, CAi, is obtained from (|2.ip - (|2.2l) - (|2.4[) and produces the discretized adjoint matrix 
(|A 2[) . The second set, CA2, is obtained from (|2.1[) - (|2.2[) - (|2.5j) and produces the discretized 
adjoint matrix (|A 3[) . The difference arises merely because the governing equations are 
arranged differently. It has no physical significance. For the DA method ( £I6.3[) the adjoint 
is simply the negative Hcrmitian of the direct matrix: <P = —T H . 

The DA method has the same truncation errors as the discretized direct system, while 
methods CAx and CA 2 have different truncation errors. The effect of these truncation 
errors is quantified in figure[Tl which compares the discrepancy between CKi and DA with 
the discrepancy between CA 2 and DA. Method GAi has generally a greater discrepancy 
than CA2, as shown in fig.[U This discrepancy is a function of the time-delay, r and the 
damping coefficients, c± and C2- Regardless of the value of the damping, the discrepancy 
is zero when r = 0. This can be inferred by examining the mathematical structure of the 
matrices, given in eq. (|A lj) - (|A 2| and (jA 3|) . If r = then <P = —T H regardless of the 
formulation used. 
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Table 1. The two different sets of continuous adjoint equations. 
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Figure 1. The discrepancy between the discrete adjoint (DA) and the continuous adjoint (CA) 
discretizations, for the two different formulations of the continuous adjoint equations, CAi and 
CA 2 . In both figures, N = 10, x h = 0.25, and /3 = 0.5. The left figure has a = 0.01, c e = 0.004 
and the right figure has eg = 0, r = 0.01. 



These adjoint equations govern the evolution of the adjoint variables, wh ich can be re- 

garde d as Lagrange multipliers from a constrained optimization perspective (jBelegundu &: Arora 
1985). Therefore, u + is the Lagrange multiplier of the acoustic momentum equation (12. ip . 
Physically, it reveals the spatial distribution of the system's sensitivity to a force. Like- 
wise, p + is the Lagrange multiplier of the pressure equation (|2.2p & (|2.4p as well as (|2.2I) 
& (|2.5p . Physically, it reveals the spatial distribution of the system's sensitivity to heat 
injection. 



5. Modal analysis: the eigenvalue problem 

So far we have considered the thermo-acoustic system in the {x, t) domain. In modal 
analysis, we consider it in the (x, a) domain using the transformations 

u(x, t) — u(x, o-)e at , u + (x, t) — ii + (x, a)e~ at , (5.1) 



p(x,t) =p{x,o-)e°\ p + (x,t)=p+(x,a)e- st . (5.2) 

where the symbol * denotes an eigenfunction. The behaviour of the system in the long 
time limit is dominated by the eigenfunction whose eigenvalue has the highest growth 
rate. The complex conjugate adjoint eigenfunctions of velocity and pressure are labelled 
u and p, respectively. With the definition of the Green's identity (|4.ip . the adjoint eigen- 
values, —a, are the negatives of the complex conjugates of the direct eigenvalues, a. 
This satisfies the bi-orthogona lity condition between the direct and adjoint eigenfunc- 
tions ( Salwen fc Grosch Il98ll) . The system is studied in the complex frequency domain 
by substituting the relations (I5.ip - (I5.2P into the direct equations (|2.ip - (l2.2p - (|2.4p and 
into the adjoint equations given in table [TJ 

Figure [2] shows the direct eigenfunctions and figure [3] the DA adjoint eigenfunctions 
as j3 increases from to 0.5. When (3 = 0, the eigenfunctions are the natural acoustic 
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Figure 2. The direct eigenfunctions as a function of the heat-release parameter, ft, for 
N = 10, x h = 0.25, r = 0.01, a = 0.01 and c 2 = 0.004. The relevant eigenvalues are: 
a = -0.0070 + 3.1416i, for j3 = 0; a = -0.0056 + 3.1848i, for /3 = 0.1; a = +0.00023 + 3.3570i, 
for /3 = 0.5. Note that the top left and bottom right frames have very small vertical scales. 
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Figure 3. The adjoint eigenfunctions as a function of the heat-release parameter, ft, for N — 10, 
Xh = 0.25, r = 0.01, ci = 0.01 and eg = 0.004. Note that the top right and bottom left frames 
have very small vertical scales. 



modes of the duct but, as j3 increases, the eigenfunctions become distorted by the heat 
release at the wire. This has important consequences for the structural sensitivity, as will 
be shown in fJT] 

Figure [4] shows the direct and adjoint eigenfunctions, found using the DA, CAi, and 
CA2 methods, at j3 = 0.5. This is the value of j3 used for the sensitivity analyses. The 
discrepancies in Im( u + ) and Re(p + ) cause the differences in sensitivities seen in £17.11 and 
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Figure 4. The adjoint eigenfunctions found using the DA, CAi, and CA2 methods. The pa- 
rameters are N = 10, x h = 0.25, r = 0.01, /3 = 0.5, cj = 0.01 and c 2 = 0.004. Note that the top 
right and bottom left frames have very small vertical scales. 



6. Calculation of the structural and base-state sensitivities 

6.1. Structural sensitivity via the CA method 
The ther mo-acoustic system describe d in $2] has been linearized about a base state. 



Following iGiannetti fc Luchini I (|2007l h we perturb the linearized operator, L, by adding 
to it some general function of the perturbation state variables, u and p. In this section, 
we assume that this feedback does not affect the base state. We also assume that the 
structural perturbation is small enough for the new thermo-acoustic configuration to be 

O-new — & + Sa, p new — p + Sp, U new = U + SH, (6.1) 

where Scr — ecr, Sp — ep, Su = eu with |e| 1, and where terms of order e 2 are sufficiently 
small to be neglected. 

The direct eigenfunctions can be arranged as a column vector [u p] T . In general, a 
structural perturbation to the thermo-acoustic operator can be represented by a 2 x 2 
tensor, SH, which acts on [u p] T . Each component SHy of this structural perturbation 
tensor quantifies the effect of a feedback mechanism between the j th eigenf unction and 
the i th direct governing equation. 

We obtain the eigenvalue drift, 5o~, caused by the structural perturbation, SH, by 
applying the Gree n's identity (14.111 to the per turbed direct and adjoint equations, in a 



manner similar to Giannetti fc Luchini | (|2007l) . Tabled describes the effect of a generic 



perturbation SH. The great advantage of this approach is that, once the direct and adjoint 
eigenfunctions have been calculated, all linear feedback mechanisms can be examined at 
little extra cost. 

We will illustrate the process for the specific case where the feedback mechanism is a 
second hot wire, called the control wire, whose parameters are denoted with the subscript 
c. The structural perturbation caused by the control wire is represented by the tensor in 
table [3l The component SH21 represents a feedback mechanism that is proportional to 
the velocity perturbation and that perturbs the pressure equation. The component SH22 
represents a feedback mechanism that is proportional to the pressure perturbation and 
that also perturbs the pressure equation. The change in the eigenvalue caused by the 
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Method 
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CA 2 




5a = 


J L [u+ p+]SH[u 


p] T dx 


J L [u+ p+]SH[u 


p] T dx 


J L (uu + +pp + )dx- 


-0TU h Ph 


J L (uu++pp+ 


)dx 



Table 2. The eigenvalue drift caused by a generic structural perturbation, which is represented 
by the generic tensor 8H. The two methods, CAi and CA2, are derived from two equivalent 
versions of the governing equations, §3] L is the dimensionless tube length. 



Method || 




CAi | 


CA 2 


«- 1 


sp B (i - 


0" 

(TT C )S(X — X c ) 








5/3 c 5(x - Xc) Sp c T c S(x - x c )-§^ 





Table 3. The tensor representing a structural perturbation caused by a second hot wire. Two 
representations are obtained, depending on whether the heat-release rate is expressed following 
the CAi or CA 2 method. 
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J L (uU + +pp + )dx+j3ruh0l 
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Table 4. The change in the eigenvalue due to the presence of the control wire with a small 
heat-release parameter 6f3 c , derived via the CAi and CA 2 approaches. 



presence of the control hot wire with a small heat-release parameter 5(3 C is given in table 
[Hfor both CA methods. The results will be described in section [7^1 

6.2. Base-state sensitivity via the CA method 

Using adjoint techniques, a single calculation can reveal how the growth rate and fre- 
quency of the system is altered by any small variation of the base-state parameters 6/3, 
(5£, St, and Sxh- This is known as the base-state sensitivity. In this section, we calcu- 
late the base-state sensitivities to j3, r, and £ as functions of the hot wire position, Xh- 
By applying a methodology similar to that presented in §6.14 we obtain the base-state 
sensitivities shown in table [5] The results will be described in §7.41 

6.3. Both sensitivities via the DA method 

Both sensitivities can be calculated directly from the discretized governing equations (the 
DA method). There are four stages to this method: (1) calculate the perturbation matrix 
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Table 5. The change in the eigenvalue due to small changes in the base-state coefficients, 
derived via the CAi and CA2 methods. 



5P using (|A4[) . imposing an arbitrarily small perturbation on the base-state parameter; 
(2) calculate the eigenvectors of the matrices f and — f H \ (3) apply (|6.2|) to find the 
eigenvalue drift; (4) divide the eigenvalue drift by the perturbation used in stage 1 in 
order to obtain the sensitivity coefficien t. The eigenvalue drift due t o a perturbation of 
the discretized direct system (similar to Giannetti fc Luchini ( 20071 )) is given by 



5a 



(6.2) 



where the matrix 5P represents a small perturbation to the direct system, whose matrix 
is T. Here, the symbol " represents an eigenvector. The column vector £ is the eigenvector 
of the adjoint matrix <P = —T H . The perturbation matrix 5P is described in (|A 4j) . It 
can represent either a structural perturbation or a base-state perturbation. 



7. Results and physical interpretation 

7.1. Comparing the three methods of calculating the structural sensitivity 

The top frames of figure [5] show the real and imaginary components of 5a/5/3 c as a 
function of the control wire position, x c , via the DA, CAi and CA2 methods. In this case 
the main hot wire i s placed at x = 0.25 so that most of the perturbation energy is in the 
first acoustic mode (Matvcev 2003). Similarly, the top frames of figure|6]show the real and 
imaginary components of 5a / 5j3 c as a function of the control wire position, x c , via the DA, 
CAi and CA 2 methods. In this case the main hot wire is pl aced at x = 0.6 25 so that most 
of the perturbation energy is in the second acoustic mode ([Matveev 1120031 ) . These results 
can be compared with the exact solution, which is obtained by finite difference. This is 
the difference between the dominant eigenvalues of the perturbed direct matrix, T + 5P, 
and the original direct matrix, T, divided by the (finite) arbitrarily small perturbation. 
The perturbation matrix 5P is given in (|A 4|) . 

As expected, the DA method matches the finite difference method exactly. The CA 
methods both contain some error, due to the truncation errors in the discretization. 
The CA2 method is usually more accurate than the CAi method. For this thermo- 
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Figure 5. Top frames: sensitivity of the growth rate, Ke(8a/Sp c ), and of the angular frequency, 
lm(8a/8/3 c ), when a control wire is placed at position x c . This is calculated exactly, via finite 
difference, and via the DA, CAi and CA2 methods. (The DA method gives the same result 
as the finite difference method to machine precision.) Bottom frames: the Rayleigh index for 
a control wire placed at x c . The parameters are N — 10, /3 — 0.5, cj = 0.01, eg = 0.004 and 
t — t c — 0.01. The main hot wire is at Xh = 0.25 so that the first acoustic mode is excited. 



acoustic system, however, the DA method turns out to be the most accurate and easy to 
implement. 

The real component of the structural sensitivity gives the change in the growth rate 
that is caused by the control wire. The imaginary component gives the change in the 
frequency. The physical reason for these changes is given in i)7.2l The control wire has a 
much stronger effect on the frequency than on the growth rate, for reasons given in i)7.3l 



7.2. Comparing the structural sensitivity with the Rayleigh Index 



It has long been known (| Rayleigh 1 118781 ) that if pressure and heat-release fluctuations 
are in phase, then acoustic vibrations are encouraged. More precisely, the Rayleigh cri- 
terion states that the energy of the acoustic field grows over one cycle of oscillation if 
§t Jt>PQ dZ^dt , exceeds the damping, where T> is the flow domain and T is the period. 
It is particularly informative to plot the spatial distribution of 

pq dt (7.1) 

T 

which is known as the Rayleigh Index. This reveals the regions of the flow that contribute 
most to the Rayleigh Criterion and therefore gives insight into the physical mechanisms 
that alter the amplitude of the oscillation. To examine the effect of the control wire, we 
substitute the approximate expressions p = pexp(ait) and q = qexp(ait) (found from 
2.4l or [2"3|) into (|7.1[) and integrate over a period 2n/<Ji, where Oi = Im(cr). (The approx- 
imation arises because the growth rate over the cycle has been ignored.) The real part of 
the Rayleigh Index gives the change in the growth rate and the imaginary part gives the 
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Figure 6. Top frames: sensitivity of the growth rate, Ke(8a/Sp c ), and of the angular frequency, 
lm(8a/8/3 c ), when a control wire is placed at position x c . This is calculated exactly, via finite 
difference, and via the DA, CAi and CA2 methods. (The DA method gives the same result 
as the finite difference method to machine precision.) Bottom frames: the Rayleigh index for 
a control wire placed at x c . The parameters are N — 10, /3 — 0.5, cj = 0.01, eg = 0.004 and 
t — t c = 0.01. The main hot wire is at xu = 0.625 so that the second acoustic mode is excited. 



change in the frequency ( figure [5]|6]). As expected, the sign of the Rayleigh index matches 
that of the structural sensitivity (the position at which it is zero matches within 1%) 
and the shape is similar. The Rayleigh Index physically explains the effect of adding the 
control hot wire to the Rijke tube. 

Firstly, we refer to fig. [5] where the main hot wire is at Xh — 0.25 and most of the 
perturbation energy is contained in the first mode. For x = to 0.56, the pressure and 
heat release eigenfunctions are sufficiently in phase that the contribution to growth over 
a cycle is positive. For x = 0.56 to 1, they are out of phase so their contribution to growth 
over a cycle is negative. For this case, the location where the presence of a second hot 
wire is most effective at reducing the growth rate is x c ~ 0.8. It is interesting to note that 
this system becomes more unstable when the control wire is placed at 0.5 < x c < 0.56. 
This is in the second half of the tube and, in the absence of the first hot wire, a control 
wire placed here would be stabilizing. The reason for this is that the main hot wire, at Xh, 
causes the eigenfunctions to distort from the acoustic modes of the duct. In particular, 
the features of the u and p eigenfunctions (figured]) shift down the duct, to higher values 
of x. This shifts downstream the region in which the control wire is destabilizing. 

Secondly, we refer to figure[B]where the main hot wire is at Xh — 0.625 and most of the per- 
turbation energy is contained in the second mode. For < x < 0.23 and 0.47 < x < 0.77, 
the pressure and heat release eigenfunctions are sufficiently in phase that the contribu- 
tion to growth over a cycle is positive; for 0.23 < x < 0.47 and 0.77 < x < 1, they are 
out of phase so their contribution to growth over a cycle is negative. For this case, the 
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location where the presence of a second hot wire is most effective at reducing the growth 
rate is x c ~ 0.36. 

7.3. Using the structural sensitivity to find the most efficient feedback mechanisms 

In passive control, an object that is placed at a point in the system causes feedback at that 
point. Under these conditions, the structural sensitivity reveals the feedback mechanism 
that is most effective at changing the frequency or growth rate of the system. 

In ij6.ll we defined the perturbation tensor, SH, to be an operator localized at the con- 
trol wire's location. In this section we consider the case of a generic feed back mechanism, 



repres ented by a localized perturbation in which SH is constant, following lGiannetti fc Luchini 

For clarity, we re-label SH as SC for this case. The structural sensitivity tensor 
S = Sa/SC is then given by the expression in table [6j Its numerator is the dyadic product 
[u + p + ] T (g> [u p] T ■ The four components of S quantify how a feedback mechanism that 
is proportional to the state variables affects the growth rate and frequency of the system. 
They are shown in fig. [7] (real part) and fig. [8] (imaginary part) as a function of x, which 
is the location of the structural perturbation. The eigenfunctions are calculated with 
both 10 modes (thick line) and 100 modes (thin line). With the latter discretization it is 
possible to capture the eigenfunction discontinuity at the hot wire's location caused by 
the impulsive heat release. Although a discretization with 100 modes does not meet the 
physical constraint that r <C 2/N O, we can use it to examine the numerical accuracy 
of the 10-mode discretization. At the hot wire's location, t he 100-mode discretization 
of Re(Si2) and Im (Su) experiences the Gibbs phenomenon ( Gibbs 18981 ) therefore the 



solution is inaccurate. The Gibbs phenomenon remains as the number of Galerkin modes 
increases. The 10-mode discretization is very accurate except at the discontinuity at the 
hot wire's location. 

When (3 — > 0, the direct eigenfunctions are nearly the acoustic modes of the system, 
as shown in fig. EE] By inspection of these eigenfunctions, we see that Su ~ (cos to) 2 , 
S12 ~ — i(sinTO) x (costo), S21 ~ i(sin7Tir) x (costo), and S22 ~ (suito) 2 , when /3 — > 0. 
The heat release from the main hot wire distorts these eigenfunctions (figure [2]) so the 
structural sensitivities are similarly distorted. 

Firstly, we consider a feedback mechanism that is proportional to the velocity and 
that forces the momentum equation {Su). For example, this could be caused by the drag 
from a fine mesh placed in the flow. The system is most sensitive when this feedback 
mechanism is placed at the entrance or exit of the duct. This is because (i) the velocity 
mode is maximal there and (ii) the adjoint velocity, which is a measure of the sensitivity 
of the momentum equation, is also maximal there, as shown in figure 01 The Re(Su) 
component (fig. [7]) is positive for all values of x, which means that, whatever value of x 
is chosen, the growth rate will decrease if the forcing is in the opposite direction to the 
velocity, as it would be for a fine mesh placed in the flow. This type of feedback greatly 
affects the growth rate (fig. [TJ) , which is the real component of the sensitivity, but barely 
affects the frequency (fig. [5]), which is the imaginary component. This behaviour is as 
expected for this type of feedback. 

Secondly, we consider a feedback mechanism that is proportional to the pressure and 
that forces the pressure equation (£22). This type of feedback is described in lChu ( 1963h 



and is relevant to pressure-coupled heat release in solid rocket engines. For this feedback, 
the system is most sensitive around the centre of the duct, with a maximum at x « 0.58. 
Again, this feedback greatly affects the growth rate (fig. [7]), and it is positive for all 
values of x, but barely affects the frequency (fig. [8]). If the heat release increases with the 
pressure, as it does for most chemical reactions, this feedback mechanism is destabilizing. 
Thirdly, we consider S12, which represents feedback from the pressure into the mo- 
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Table 6. Structural sensitivity tensor for a general feedback mechanism SC. 
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Figure 7. Real part of the components of the structural sensitivity tensor (via the DA method) 
for the same parameters as fig. [5] These components indicate the effect of a feedback mechanism 
on the growth rate of oscillation. 



mentum equation and S21 , which represents feedback from the velocity into the pressure 
equation. These types of feedback barely affect the growth rate (fig. [7j) but greatly affect 
the frequency (fig. [8]). The hot control wire considered in figure [5] causes this type of 
feedback (S21) if t — 0. This analysis shows, therefore, that this passive control device 
is quite ineffective at reducing the growth rate. This had been seen already in figure 
in which the hot wire is seen to affect the frequency (imaginary component) much more 
than it affects the growth rate (real component). 

By inspection of fig. [71 we conclude that the passive device that is most effective at 
reducing the growth rate should force the momentum equation in the opposite direction 
to the velocity fluctuation and should be placed at the exit of the tube. A damping device 
such as an adiabatic wire mesh would achieve this. 

This paper is mainly relevant to passive control but it is worth briefly mentioning 
active control. For active control, the sensor and actuator would typically be in different 
places. For maximum observability, the sensor should be placed where the relevant direct 
eigenfunction has its largest amplitude. For maximum controllability, the actuator should 
be placed where the relevant adjoint eigenfunction has its largest amplitude. 
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Figure 8. Imaginary part of the components of the structural sensitivity tensor (via the DA 
method) for the same parameters as fig. [5] These components indicate the effect of a feedback 
mechanism on the angular frequency of oscillation. 



7.4. Base-state sensitivity results 

Figure^ shows how a small variation in the heat-release parameter, /3, affects the growth 
rate, Re(er), and the angular frequency, Im(cr), for different hot-wire positions, x^. Figure 
[Hb shows how a small variation in the time-delay coefficient, r, affects the same quantities. 
These are calculated via the DA, CAi, and CA2 methods and the results are checked 
against the exact solution, which is obtained by finite difference, as in £16.11 

As shown in tabled these curves depend on the shapes of the direct and adjoint eigen- 
functions. In turn, these eigenfunctions are distorted from the natural acoustic modes of 
the duct by the heat release from the wire. (This distortion is shown in figures [2] and [3] 
for xt = 0.25.) This accounts for the elaborate shapes of the base flow sensitivity curves. 
It is also worth commenting on their relative magnitudes: small variations in (3 have a 
much greater effect on the frequency than on the growth rate, while small variations in r 
have a much greater effect on the growth rate than on the frequency. This will always be 
the case when lot <C 1, which is easy to justify by the following argument. If p ~ sinwt at 
the hot wire, then u ~ cos tut and q ~ coscj(£ — r) there. Using trigonometric relations, 
it is easy to show that § pq dt, which quantifies how much ft affects the growth rate, 
is proportional to sinwr and that § uq dt, which quantifies how much j3 affects the fre- 
quency, is proportional to coswr. Therefore, for small lot, the change in the growth rate, 
Re(5a/5/3), should be of order lot, while the change in the frequency, lm(Sa/S/3), should 
be of order 1. Differentiating with respect to r at constant /3, we find that the change in 
f pq dt due to a change in r is proportional to lo cos lot. Similarly, the change in § uq dt 
due to a change in r is proportional to losycilot. Therefore, for small lot, Re(<5(7/<5r) 
should be of order lo, while lm(Sa/Sf3) should be of order lo 2 t. These magnitudes closely 
match the amplitudes in figure |H1 for which lo « 7t and r = 0.01. 

Figure |9p shows how the angular frequency changes with the damping factor £. A small 
increase in £ lowers the frequency of the linear oscillations. A small increase of £ is always 
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stabilizing, i.e. the growth rate decreases, but does not depend on the hot-wire position 
(figure not shown). In order to study the sensitivity to small changes of the damping, 
<5£, only one Galcrkin mode has been considered. This is because £ is a function of the 
Galerkin mode, as explained in fJ3] Therefore, with the damping model and numerical 
discretization adopted, formulae in the bottom row of table [5] are valid only for the first 
Galerkin mode. 

As for the structural sensitivity, there is a discrepancy between the DA and CA solu- 
tions, which arises from the different truncation errors in the discretizations. The origin 
of this error can be inferred from the matrices in appendix A. The CAi method provides 
an inaccurate lva{8a / 5t) 1 as shown in figure[9b. This is due to the time-delay coefficient 
and this discrepancy vanishes as the time-delay becomes much smaller. In this case, we 
find that the maximal discrepancy between CAi and the exact solution is smaller than 
10% when t < 0.001. 



8. Passive control of an unstable system 

In this section we demonstrate the suppression of thermo-acoustic oscillations using 
a control wire placed at the optimal location, as predicted by the structural sensitivity 
analysis. We use the parameters in figure [5l which shows that, in order to reduce the 
growth rate most effectively, the control wire should be placed at x c = 0.8. We integrate 




Figure 10. Stabilization ol the thermo-acoustic system via a second hot wire introduced at 
t = 1000 and x c — 0.8. f3 c = /3/10 = 0.05 and the remaining parameters are the same as in fig. [5] 
The time integration (a) is performed on the nonlinear time-delayed equations discretized with 
20 Galerkin modes. The solution is shown at x = 0.25. The amplitude of the spectrum of the 
solution is shown in (b)-(c). 



the nonlinear time-delayed governing equations (|B ![) - (|F3 2j) forward in time with a 4 th 
order Runge-Kutta algorithm and 20 Galerkin modes. 

When the control wire is absent, the growth rate is ay = 0.00023 and the angular 
frequency is <7j = 3.3570. We set the heat-release parameter for the control wire to be j3 c = 
/3/10 = 0.05, which is small enough to fulfil the linear assumptions. When the control wire 
is present, the growth rate is ay = —0.00058 and the angular frequency is <7j = 3.3354. 
The difference between these values matches that predicted by the structural sensitivity 
analysis, for which da = p c xSa/S/3 c 0.05 x (-0.01633-0.4323i) = -0.00082-0.02162i, 
at x c — 0.8. 

Figure [TTJa shows the pressure at x = 0.25 as a function of time in the nonlinear 
simulations. The control wire is introduced at t = 1000. The behaviour is as expected: 
there is exponential growth until t = 1000 and exponential decay afterwards. In fig. 1 10b - 
110b the fast Fourier transform (FFT) performed on the nonlinear time-solution confirms 
the frequency shift predicted by the sensitivity analysis. 



9. Conclusions 

The main goal of this paper is to take a technique developed for the analysis of hydro- 
dynamic stability and adapt it to the analysis of thermo-acoustic stability. This technique 
uses adjoint equations to calculate a system's sensitivity to feedback or to changes in the 
base state. 

By arranging the linearized thermo-acoustic governing equations in two different ways, 
we derive two different sets of adjoint equations, which we then discretize with a Galerkin 
decomposition. This is known as the 'Continuous Adjoint' (CA) method and the two 
sets of adjoint equations produce two different matrices, labelled CAi and CA 2 . We 
also derive the adjoint equations directly from the discretized linearized thermo-acoustic 
system. This is known as the 'Discrete Adjoint' (DA) method and it produces another 
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matrix, labelled DA. The DA matrix is the negative Hcrmitian of the matrix representing 
the discretized governing equations. We calculate the direct and adjoint eigenfunctions 
of the thermo-acoustic system using these direct and adjoint matrices. We find that the 
DA method is more accurate and easier to implement than either CA method for this 
thermo-acoustic model. 

Two sensitivity analyses are carried out: one focuses on structural perturbations and 
the other on base-state perturbations. In the structural sensitivity analysis, we calculate 
the effect that a generic feedback mechanism has on the frequency and growth rate of 
oscillations. We illustrate this by considering the influence of a second hot wire, with a 
small heat release parameter. We find that the second wire affects the frequency much 
more than the growth rate and explain this physically by evaluating the Rayleigh Index 
for the second hot wire. We then use the results of the structural sensitivity to identify the 
feedback mechanism that is most effective at reducing the growth rate of oscillations. We 
find that this mechanism should force the momentum equation in the opposite direction 
to the velocity perturbation and that it should be placed at the downstream end of the 
duct. An adiabatic fine mesh would achieve this. In the base-state sensitivity analysis, we 
calculate the effect that a small variation in the base-flow parameters has on the frequency 
and growth rate of oscillations. As expected, we find that a small increase in the wire 
temperature affects the frequency more than the growth rate and that a small increase in 
the time delay affects the growth rate more than the frequency. Also as expected, we find 
that a small increase in the damping always has a stabilizing effect. The novelty of this 
paper is in the technique. Each sensitivity analysis was obtained extremely quickly with 
a single calculation. It was then checked against the exact solution found by many finite 
difference calculations. The DA method matched the finite difference method exactly, 
while there was some discrepancy when using the CAi and CA2 methods. 

The successful application of sensitivity analysis to thermo-acoustics opens up new 
possibilities for the passive control of thermo-acoustic oscillations. In a single calcula- 
tion, sensitivity analysis shows how the growth rate and frequency of small oscillations 
about some base state are affected either by a passive control element in the system or 
by a change to its base state. This gradient information can be combined with other con- 
straints, such as that the total mean heat release be constant, to show how an unstable 
thermo-acoustic system should be changed in order to make it stable. In this paper, we 
have demonstrated this for a simple system with four elements to the base state: the hot 
wire position, its heat-release coefficient, its time delay and the damping. In future work, 
we will examine more elaborate flame models and acoustic networks. This will allow us 
to calculate the sensitivity to the flame shape and to the characteristics of the acoustic 
network in which the flame sits. 

We would like to thank Dr. Outi Tammisola (University of Cambridge, Department of 
Engineering, U.K.) for valuable discussions and comments on this paper. This work was 
supported by the European Research Council through Project ALORS 2590620. 

Appendix A. Discretized equations 

It is useful to define the following matrices, which are expressed in matrix notation 
(repeated indices are not to be summed): 

Aj = 0, B,j = nSiji, Eij(a,c 2 ) = -C&j, 
Fij{Pw,x w ) = —2j3 w sin(TOXtu) cos(njx w ), Gij(/3 w ,x w ,T w ) = 2inT w (3 w sm(nix w ) cos(njx w ), 
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) = 2/3 w r w (j cos(nix w ) sm(njx w ), 

Q(^1l?l &w) = By F Fjj , Djj({3 w , X W7 T w , Cl, C2) = ~t~ Gy. 

where i, j = 1,2, ...,N, N is the number of Galerkin modes, <5y is the Kronecker delta 
and w stands for wire. The direct matrix V is given by: 



r 



>4 e 

C(/3 h ,x h ) D(j3 h ,x hl T h ,ci,c 2 ) 



(Al) 



The continuous adjoint equations (table [IJ are discretized using the Galerkin method 
as for the direct modes, by means of the decomposition in eq. (|3.1[1 . The discretization of 
the first set of adjoint equation CAi (table [1]) gives rise to the following adjoint matrix 







-G T (/3 hl x h ,T h ) 
-B 



B - F T (f3 hl x h ) + H(P h ,x h ,Th,ci,c 2 ) 



-E(ci,c 2 ) 

while the second set CA2 (table [IJ gives the following adjoint matrix 



B 



B-F T ((3 h ,x h ) 
-£(ci,c 2 ) + G T (/3 h ,x h ,T h ) 



(A 2) 



(A3) 



Note that T and are 2N x 2N matrices. We indicated the main hot wire with subscript 
h and the control hot wire with the subscript c. Finally, the perturbation matrix of the 
direct system is: 



5P = 



[0] 



NxN 



[0] 



NxJV 



Sf3 h ,x h ) + . 
- C(6f3 c ,x c ) 



D(/3 h + S/3 h ,x h ,T h +St, 1 ,c 1 +Scx,c 2 
... + D(8i3 c ,x c ,5t c ,ci,c 2 ) 



5c 2 ) 



(A4) 

On the one hand, we obtain the perturbation matrix caused by the presence of the second 
hot wire by setting <5/3/ t = St^ = 5cj = Sc2 =0 and 8f5 c > and 6t c > 0. On the other 
hand, we obtain the perturbation matrix caused by (positive) base-flow variations by 
setting 5f3 c = St c = and 5(3h > 0, Stu > 0, Scj > and Sc2 > 0. 



Appendix B. Nonlinear time- delayed equations for control 

In this section we provide the nonlinear time-delayed equations of the thermo-acoustic 
system with a control hot wire. 

Referring to the time integration presented in when the second hot wire is off, for 
t < 1000, then (3 C = 0; when the second wire is on, for t > 1000, then (3 C = (3/10. 



du dp 
dt dx 



(Bl) 




Xh, 



(B2) 
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